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Abstract 

We build up the mathematical connection between the "Expectation-Maximization" (EM) algorithm and 
gradient-based approaches for maximum likelihood learning of finite Gaussian mixtures. We show that the 
EM step in parameter space is obtained from the gradient via a projection matrix P , and we provide an 
explicit expression for the matrix. We then analyze the convergence of EM in terms of special properties 
of P and provide new results analyzing the effect that P has on the likelihood surface. Based on these 
mathematical results, we present a comparative discussion of the advantages and disadvantages of EM 
and other algorithms for the learning of Gaussian mixture models. 
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1 Introduction 

The "Expectation-Maximization" (EM) algorithm is a 
general technique for maximum likelihood (ML) or max- 
imum a posteriori (MAP) estimation. The recent em- 
phasis in the neural network literature on probabilistic 
models has led to increased interest in EM as a possible 
alternative to gradient-based methods for optimization. 
EM has been used for variations on the traditional theme 
of Gaussian mixture modeling (Ghahramani k Jordan, 
1994; Nowlan, 1991; Xu k Jordan, 1993a, b; Tresp, Ah- 
mad k Neuneier, 1994; Xu, Jordan k Hinton, 1994) and 
has also been used for novel chain-structured and tree- 
structured architectures (Bengio k Frasconi, 1995; Jor- 
dan k Jacobs, 1994). The empirical results reported in 
these papers suggest that EM has considerable promise 
as an optimization method for such architectures. More- 
over, new theoretical results have been obtained that 
link EM to other topics in learning theory (Amari, 1994; 
Jordan k Xu, 1993; Neal k Hinton, 1993; Xu k Jordan, 
1993c; Yuille, Stolorz k Utans, 1994). 

Despite these developments, there are grounds for 
caution about the promise of the EM algorithm. One 
reason for caution comes from consideration of theoret- 
ical convergence rates, which show that EM is a first 
order algorithm. 1 More precisely, there are two key re- 
sults available in the statistical literature on the con- 
vergence of EM. First, it has been established that un- 
der mild conditions EM is guaranteed to converge to 
a local maximum of the log likelihood / (Boyles, 1983; 
Dempster, Laird k Rubin, 1977; Redner k Walker, 
1984; Wu, 1983). (Indeed the convergence is monotonic: 
/(0( & + 1 )) > 1{Q^)), where 6^ is the value of the pa- 
rameter vector at iteration k.) Second, considering 
EM as a mapping 6 (&+1) = M(6 (&) ) with fixed point 
6* = M(0*), we have e( &+1 )-6* « dM d ( ®,' ] (& (k) -&*) 
when Q( k+1 > is near 0*, and thus 
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almost surely. That is, EM is a first order algorithm. 

The first-order convergence of EM has been cited in 
the statistical literature as a major drawback. Red- 
ner and Walker (1984), in a widely-cited article, argued 
that superlinear (quasi-Newton, method of scoring) and 
second-order (Newton) methods should generally be pre- 
ferred to EM. They reported empirical results demon- 
strating the slow convergence of EM on a Gaussian mix- 
ture model problem for which the mixture components 
were not well separated. These results did not include 
tests of competing algorithms, however. Moreover, even 
though the convergence toward the "optimal" parameter 
values was slow in these experiments, the convergence in 
likelihood was rapid. Indeed, Redner and Walker ac- 
knowledge that their results show that "... even when 



An iterative algorithm is said to have a local convergence 
rate of order q > 1 if ||0(* +1 ) - Q*\\/\\Q^ - 0*f < r + 
o(\\& (k) - 0*||) for k sufficiently large. 



the component populations in a mixture are poorly sep- 
arated, the EM algorithm can be expected to produce 
in a very small number of iterations parameter values 
such that the mixture density determined by them re- 
flects the sample data very well." In the context of the 
current literature on learning, in which the predictive 
aspect of data modeling is emphasized at the expense of 
the traditional Fisherian statistician's concern over the 
"true" values of parameters, such rapid convergence in 
likelihood is a major desideratum of a learning algorithm 
and undercuts the critique of EM as a "slow" algorithm. 

In the current paper, we provide a comparative anal- 
ysis of EM and other optimization methods. We empha- 
size the comparison between EM and other first-order 
methods (gradient ascent, conjugate gradient methods), 
because these have tended to be the methods of choice 
in the neural network literature. However, we also com- 
pare EM to superlinear and second-order methods. We 
argue that EM has a number of advantages, including its 
naturalness at handling the probabilistic constraints of 
mixture problems and its guarantees of convergence. We 
also provide new results suggesting that under appropri- 
ate conditions EM may in fact approximate a superlin- 
ear method; this would explain some of the promising 
empirical results that have been obtained (Jordan k Ja- 
cobs, 1994), and would further temper the critique of EM 
offered by Redner and Walker. The analysis in the cur- 
rent paper focuses on unsupervised learning; for related 
results in the supervised learning domain see Jordan and 
Xu (in press). 

The remainder of the paper is organized as follows. 
We first briefly review the EM algorithm for Gaussian 
mixtures. The second section establishes a connection 
between EM and the gradient of the log likelihood. We 
then present a comparative discussion of the advantages 
and disadvantages of various optimization algorithms in 
the Gaussian mixture setting. We then present empir- 
ical results suggesting that EM regularizes the condi- 
tion number of the effective Hessian. The fourth section 
presents a theoretical analysis of this empirical finding. 
The final section presents our conclusions. 

2 The EM algorithm for Gaussian 
mixtures 

We study the following probabilistic model: 
K 



p(x\e) = J2»sP(*\™s^ j ), 



(i) 



P{x\m h T,j) 



(2tt)' j / 2 |Sj| 1 / 2 



■t;(x — mj) S. (x — nij) 



where oij > and X2j=i a i = 1 anc ^ d ls ^ ne dimension- 
ality of the vector x. The parameter vector consists 
of the mixing proportions oij , the mean vectors rrij , and 
the covariance matrices T,j . 

Given K and given N independent, identically dis- 
tributed samples {x^'}^ , we obtain the following log 



likelihood: 
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/(e) = io g n^ (t) |e) = ^iogP(^)|e), (2) 
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which can be optimized via the following iterative algo- 
rithm (see, e.g, Dempster, Laird & Rubin, 1977): 
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where the posterior probabilities h- are defined as fol- 
lows: 
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3 Connection between EM and 
gradient ascent 

In the following theorem we establish a relationship be- 
tween the gradient of the log likelihood and the step in 
parameter space taken by the EM algorithm. In par- 
ticular we show that the EM step can be obtained by 
premultiplying the gradient by a positive definite ma- 
trix. We provide an explicit expression for the matrix. 

Theorem 1 At each iteration of the EM algorithm Eq. 

(3), we have 
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where A denotes the vector of mixing proportions 
[a i, ■ ■ ■ , cykV , j indexes the mixture components (j = 
1,- ■ ■ ,K), k denotes the iteration number, "vec[B]" de- 
notes the vector obtained by stacking the column vectors 



2 Although we focus on maximum likelihood (ML) estima- 
tion in this paper, it is straightforward to apply our results 
to maximum a posteriori (MAP) estimation by multiplying 
the likelihood by a prior. 



of the matrix B, and "®" denotes the Kronecker prod- 
uct. Moreover, given the constraints E)=i a ; = 1 an d 
a- ' > 0, Pjj is a positive definite matrix and the ma- 
trices Pm- and -Pj; . are positive definite with probability 
one for N sufficiently large. 

Proof. (1) We begin by considering the EM update 
for the mixing proportions a.{. From Eqs. (1) and (2), 
we have 

N 
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A=A( k ) - 2-^i ' 
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Premultiplying by Pjj , we obtain 

0) ^ 

P A -g^\A=A( k ) 



l_ " {[a f) P{x (r K tf) )} . . .]T _ A ( k) z k =i tt (*) p(ar ( 0) gW )} 

M 2-^ 



N 



£f =ia f)p(,(*Mf)) 



4 = 1 

N 



N 



The update formula for A in Eq. (3) can be rewritten as 

1 N 
.4O+1) = A(k) + j_ £[/,(*)(*), . . . , h p(t)f _ A ( k ), 

4 = 1 

Combining the last two equations establishes the update 
rule for A (Eq. 4). Furthermore, for an arbitrary vec- 
tor u, we have Nu T P^ u = M T diag[a^ , • • -,a K ]u — 
(u T A^') 2 . By Jensen's inequality we have 
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= (u T A^) 2 . 

Thus, u T P^ u > and P^ is positive definite given 

the constraints E?=i a j 1 anc ^ Qu > for all j . 

(2) We now consider the EM update for the means 
ra,-. It follows from Eqs. (1) and (2) that 
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Premultiplying by Pm/ yields 
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From Eq. (3), we have Et=i ^j (0 ^ ^; moreover, Sy 
is positive definite with probability one assuming that N 



is large enough such that the matrix is of full rank. Thus, 

it follows from Eq. (8) that Pm] is positive definite with 
probability one. 

(3) Finally, we prove the third part of the theorem. 
It follows from Eqs. (1) and (2) that 
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With this in mind, we rewrite the EM update formula 

for £„■ ' as 
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That is, we have 
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Utilizing the identity vec[ABC] = (C T (g) A)vec[B], we 
obtain 
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Thus P™ = „ 2 , (Sf } ® Sf } ). Moreover, for an 
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arbitrary matrix U , we have 

vec[t/] T (sf ) ®sf ) )vec[t/] 

= tr(sfVEfV T ) 

= tr((sf ) P) T (EfV)) 

= vecpfV] T vecpfV] 
> 0, 

(k) 

where equality holds only when £> U = for all [/". 

Ik) 
Equality is impossible, however, since £> is positive 

definite with probability one N is sufficiently large. Thus 
it follows from Eq. (9) and ^=1 h f \ t ) > ° that p ^ 
is positive definite with probability one. □ 
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Using the notation 

6 = [m( , • • • , m T K , vec[Si] T , • • • , vecpx] 7 , .4 

and P(6) = diag[_P mi , • • • , P mK , P Sl , 

can combine the three updates in Theorem 1 into a single 

equation: 

e ^ + i) = e w + P(e w ) ^L le=e(k)j (10) 

Under the conditions of Theorem 1, P(0(*)) is a positive 
definite matrix with probability one. Recalling that for 

a positive definite matrix B, we have J^ PJfk ^ ^> we 
have the following corollary: 

Corollary 1 For each iteration of the EM algorithm 
given by Eq.(3), the search direction 0^+ 1 ^ — 0W has 
a positive projection on the gradient of I. 

That is, the EM algorithm can be viewed as a variable 
metric gradient ascent algorithm for which the projection 
matrix _P(0 ( - ') changes at each iteration as a function 
of the current parameter value Q^K 

Our results extend earlier results due to Baum and 
Sell (1968). Baum and Sell studied recursive equations 
of the following form: 
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1, where J is a polyno- 
mial in x\^' having positive coefficients. They showed 
that the search direction of this recursive formula, i.e., 
T(x^>) — x^ \ has a positive projection on the gradient 
of of J with respect to the x^ k > (see also Levinson, Ra- 
biner & Sondhi, 1983). It can be shown that Baum and 
Sell's recursive formula implies the EM update formula 
for A in a Gaussian mixture. Thus, the first statement 
in Theorem 1 is a special case of Baum and Sell's earlier 
work. However, Baum and Sell's theorem is an existence 
theorem and does not provide an explicit expression for 
the matrix Pa that transforms the gradient direction 
into the EM direction. Our theorem provides such an 
explicit form for Pj^. Moreover, we generalize Baum and 
Sell's results to handle the updates for rrij and £ j , and 
we provide explicit expressions for the positive definite 
transformation matrices P m and Py, as well. 

It is also worthwhile to compare the EM algorithm to 
other gradient-based optimization methods. Newton's 
method is obtained by premultiplying the gradient by 
the inverse of the Hessian of the log likelihood: 



0(*+i) = 0(*) + fr(0(*))- 
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Newton's method is the method of choice when it can 
be applied, but the algorithm is often difficult to use 
in practice. In particular, the algorithm can diverge 
when the Hessian becomes nearly singular; moreover, 
the computational costs of computing the inverse Hes- 
sian at each step can be considerable. An alternative 



is to approximate the inverse by a recursively updated 
matrix 5( &+1 ) = B^ + r]AB^ k \ Such a modification 
is called a quasi- Newton method. Conventional quasi- 
Newton methods are unconstrained optimization meth- 
ods, however, and must be modified in order to be used 
in the mixture setting (where there are probabilistic con- 
straints on the parameters). In addition, quasi-Newton 
methods generally require that a one-dimensional search 
be performed at each iteration in order to guarantee con- 
vergence. The EM algorithm can be viewed as a special 
form of quasi-Newton method in which the projection 
matrix P(O^) in Eq. (10) plays the role of B^ . As 
we discuss in the remainder of the paper, this partic- 
ular matrix has a number of favorable properties that 
make EM particularly attractive for optimization in the 
mixture setting. 

4 Constrained optimization and general 
convergence 

An important property of the matrix P is that the EM 
step in parameter space automatically satisfies the prob- 
abilistic constraints of the mixture model in Eq. (1). 
The domain of contains two regions that embody the 

probabilistic constraints: T>\ = {0 : X2?=i a j = 1} anc ^ 
V 2 = {0 : af ] > 0, Ej positive definite}. For the EM 
algorithm the update for the mixing proportions aj can 
be rewritten as follows: 



S k + r ) 



1 * afp(x^\nf\^) 



It is obvious that the iteration stays within V\ . Simi- 
larly, the update for T,j can be rewritten as: 
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which stays within X>2 for N sufficiently large. 

Whereas EM automatically satisfies the probabilistic 
constraints of a mixture model, other optimization tech- 
niques generally require modification to satisfy the con- 
straints. One approach is to modify each iterative step 
to keep the parameters within the constrained domain. 
A number of such techniques have been developed, in- 
cluding feasible direction methods, active sets, gradient 
projection, reduced-gradient, and linearly-constrained 
quasi-Newton. These constrained methods all incur ex- 
tra computational costs to check and maintain the con- 
straints and, moreover, the theoretical convergence rates 
for such constrained algorithms need not be the same as 
that for the corresponding unconstrained algorithms. A 
second approach is to transform the constrained opti- 
mization problem into an unconstrained problem before 
using the unconstrained method. This can be accom- 
plished via penalty and barrier functions, Lagrangian 
terms, or re-parameterization. Once again, the extra al- 
gorithmic machinery renders simple comparisons based 



on unconstrained convergence rates problematic. More- 
over, it is not easy to meet the constraints on the covari- 
ance matrices in the mixture using such techniques. 

A second appealing property of P(Q^ k >) is that each 
iteration of EM is guaranteed to increase the likelihood 
(i.e., /(©( &+1 )) > l(Q( k >)). This monotonic convergence 
of the likelihood is achieved without step-size parameters 
or line searches. Other gradient-based optimization tech- 
niques, including gradient descent, quasi-Newton, and 
Newton's method, do not provide such a simple theo- 
retical guarantee, even assuming that the constrained 
problem has been transformed into an unconstrained 
one. For gradient ascent, the step size r\ must be chosen 
to ensure that ||©( &+1 ) - e^" 1 ) | |/||(e( fc ) - O^" 1 )))) < 
||/-|-?7-ff(© ( - &-1 - ) ))|| < 1. This requires a one- dimensional 
line search or an optimization of r\ at each iteration, 
which requires extra computation which can slow down 
the convergence. An alternative is to fix ry to a very 
small value which generally makes ||7 + ry_ff(0^ _1 ^))|| 
close to one and results in slow convergence. For New- 
ton's method, the iterative process is usually required 
to be near a solution, otherwise the Hessian may be in- 
definite and the iteration may not converge. Levenberg- 
Marquardt methods handle the indefinite Hessian ma- 
trix problem; however, a one- dimensional optimization 
or other form of search is required for a suitable scalar 
to be added to the diagonal elements of Hessian. Fisher 
scoring methods can also handle the indefinite Hessian 
matrix problem, but for non-quadratic nonlinear opti- 
mization Fisher scoring requires a stepsize r\ that obeys 
\\I + rjBH(e^-^))\\ < 1, where B is the Fisher infor- 
mation matrix. Thus, problems similar to those of gra- 
dient ascent arise here as well. Finally, for the quasi- 
Newton methods or conjugate gradient methods, a one- 
dimensional line search is required at each iteration. In 
summary, all of these gradient-based methods incur ex- 
tra computational costs at each iteration, rendering sim- 
ple comparisons based on local convergence rates unre- 
liable. 

For large scale problems, algorithms that change the 
parameters immediately after each data point ("on-line 
algorithms") are often significantly faster in practice 
than batch algorithms. The popularity of gradient de- 
scent algorithms for neural networks is in part to the 
ease of obtaining on-line variants of gradient descent. 
It is worth noting that on-line variants of the EM algo- 
rithm can be derived (Neal & Hinton, 1993, Titterington, 
1984), and this is a further factor that weighs in favor 
of EM as compared to conjugate gradient and Newton 
methods. 

5 Convergence rate comparisons 

In this section, we provide a comparative theoretical dis- 
cussion of the convergence rates of constrained gradient 
ascent and EM. 

For gradient ascent a local convergence result can by 
obtained by Taylor expanding the log likelihood around 
the maximum likelihood estimate O* . For sufficiently 
large k we have: 

lie( fc + 1 )-o*||< ||/ + ?yj ff(e*))||||(o( J; )-o*)|| (12) 



|| J + T}H(e*)\\ <\ M [I+ vH(Q*)} = r, (13) 

where H is the Hessian of /, r\ is the step size, and 
r = max{|l - ^[-5(0')]!, |1 - V X m [-H(0*)}\}, 
where Xm [A] and X m [A] denote the largest and small- 
est eigenvalues of A, respectively. 

Smaller values of r correspond to faster convergence 
rates. To guarantee convergence, we require r < 1 or 
< r\ < 2/Xm[—H(Q*)]. The minimum possible value 
of r is obtained when r\ = 1/Xm[H(Q*)] with 

r-min = i-x m [H(e*)]/x M [H(e*)] 
= i-k- 1 ^*)], 

where k[H] = \M[H]/\ m [H] is the condition number of 
H . Larger values of the condition number correspond to 
slower convergence. When k[H] = 1 we have r msn = 0, 
which corresponds to a superlinear rate of convergence. 
Indeed, Newton's method can be viewed as a method 
for obtaining a more desirable condition number — the 
inverse Hessian H~ l balances the Hessian H such that 
the resulting condition number is one. Effectively, New- 
ton can be regarded as gradient ascent on a new func- 
tion with an effective Hessian that is the identity matrix: 
H e ff = H~ l H = I . In practice, however, k[H] is usually 
quite large. The larger k[H] is, the more difficult it is to 
compute H~ l accurately. Hence it is difficult to balance 
the Hessian as desired. In addition, as we mentioned 
in the previous section, the Hessian varies from point 
to point in the parameter space, and at each iteration 
we need recompute the inverse Hessian. Quasi-Newton 
methods approximate H(Q^')~ l by a positive matrix 
B^> that is easy to compute. 

The discussion thus far has treated unconstrained op- 
timization. In order to compare gradient ascent with 
the EM algorithm on the constrained mixture estima- 
tion problem, we consider a gradient projection method: 

e(H1) = eW + ^3el) < 14 ) 

where n^ is the projection matrix that projects the gra- 
dient a& (k) into T>\. This gradient projection iteration 
will remain in V\ as long as the initial parameter vector 
is in V\. To keep the iteration within X>2, we choose an 
initial 0(°) £ X>2 and keep r\ sufficiently small at each 
iteration. 

Suppose that E = [ei, • • • , e m ] are a set of indepen- 
dent unit basis vectors that span the space V\ . In this 
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ae( k ) c ~ aew ■ 
||q(*0 — ©* 1 1 . In this representation the projective gradi- 
ent algorithm Eq. (14) becomes simple gradient ascent: 
6^ +1) = 6^ } + ryT^sy. Moreover, Eq. (12) becomes 
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He^ 1 ) - e*|| < \\e t (i + r ?J ff(e*))||||(e( & ) - e*)||. As 

a result, the convergence rate is bounded by 
r c = \\E T (I + 11 H(e*))\\ 



< X X M [ET(I + 11 H{Q*)){I + 11 H{Q*))TE] 



X M [E T (I + 2r)H(e*) + r] 2 H 2 (0*))E}. 
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In this equation H c = E T H(Q)E is the Hessian of / 
restricted to V\. 

We see from this derivation that the convergence 
speed depends on k[H c ] = A«-[- H c ]/X m [— H c ]. When 
k[H c ] = 1, we have ^/l + T] 2 X 2 M (-H C ) - 2r]X m [-H c ] = 
1 — r)X[—H c ], which in principle can be made to equal 
zero if r\ is selected appropriately. In this case, a super- 
linear rate is obtained. Generally, however, k[H c ] ^ 1, 
with smaller values of k[H c ] corresponding to faster con- 
vergence. 

We now turn to an analysis of the EM algorithm. As 
we have seen EM keeps the parameter vector within V\ 
automatically. Thus, in the new basis the connection 
between EM and gradient ascent (cf. Eq. (10)) becomes 

and we have 

ye^+^-ei^ \\e t (i + PH(e*))\\\\(e ( ^ - e*)\\ 

with 

r c = \\E T (I + PH(e*))\\ 



< ^X M [E T (I + PH(Q*))(I + PH(Q*)) T E]. 
The latter equation can be further manipulated to yield: 



<\/l 



\* M [E?PHE\- 



2X m [-E T PHE]. (16) 



Thus we see that the convergence speed of EM de- 
pends on k[E t PHE] = X M [E T PHE]/X m [E T PHE]. 
When kIE^PHE] = 1, X M [E T PHE] = 1, we 
have y/l + X 2 M [E T PHE] - 2X m [-E T PHE] = (1 - 
^m[— E PHE\) = 0. In this case, a superlinear rate 
is obtained. We discuss the possibility of obtaining su- 
perlinear convergence with EM in more detail below. 

These results show that the convergence of gradient 
ascent and EM both depend on the shape of the log likeli- 
hood as measured by the condition number. When k[H] 
is near one, the configuration is quite regular, and the 
update direction points directly to the solution yielding 
fast convergence. When k[H] is very large, the / sur- 
face has an elongated shape, and the search along the 
update direction is a zigzag path, making convergence 
very slow. The key idea of Newton and quasi-Newton 
methods is to reshape the surface. The nearer it is to a 
ball shape (Newton's method achieves this shape in the 
ideal case), the better the convergence. Quasi-Newton 
methods aim to achieve an effective Hessian whose con- 
dition number is as close as possible to one. Interest- 
ingly, the results that we now present suggest that the 
projection matrix P for the EM algorithm also serves 
to effectively reshape the likelihood yielding an effective 
condition number that tends to one. We first present 
empirical results that support this suggestion and then 
present a theoretical analysis. 
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Figure 1: Experimental results for the estimation of the 
parameters of a two-component Gaussian mixture, (a) 
The condition numbers as a function of the iteration 
number, (b) A zoomed version of (a) after discarding 
the first 25 iterations. The terminology 'original, con- 
strained, and EM-equivalent Hessians' refers to the ma- 
trices H,E T HE, and E T PHE respectively. 



We sampled 1000 points from a simple finite mixture 
model given by 




p(x) = aipi(x) + a 2 p2(x) 



where 



Pi(x) 



27T<7, 



l(x-mi) 2 
2 GX P{~9 „? >• 



The parameter values were as follows: a.\ = 
0.7170, a 2 = 0.2830, m 1 = -2, m 2 = 2, a\ = 1, a\ = 
1. We ran both the EM algorithm and gradient ascent 
on the data. At each step of the simulation, we calcu- 
lated the condition number of the Hessian (k[H (0( k >)]) , 
the condition number determining the rate of conver- 
gence of the gradient algorithm (n[E T H(Q^')E\), and 
the condition number determining the rate of conver- 
gence of EM {1^ P{Q^)H{Q^)E}). We also calcu- 
lated the largest eigenvalues of the matrices H(Q^>), 
E T H(QW)E, and E T P '(6^) H (6^)) E . The results 
are shown in Fig. 1. As can be seen in Fig. 1(a), the con- 
dition numbers change rapidly in the vicinity of the 25th 
iteration and the corresponding Hessian matrices be- 
come indefinite. Afterward, the Hessians quickly become 
definite and the condition numbers converge. 3 As shown 
in Fig. 1(b), the condition numbers converge toward the 
values n[H(eW)] = 47.5, k[E t H(Q^)E] = 33.5, and 
k[E t P{Q^y)H{Q^y)E] = 3.6. That is, the matrix P 
has greatly reduced the condition number, by factors of 
9 and 15. This significantly improves the shape of / and 
speeds up the convergence. 

We ran a second experiment in which the means of the 
component Gaussians were m\ = — 1 and rri2 = 1. The 
results are similar to those shown in Fig. 1. Since the 
distance between two distributions is reduced into half, 



3 Interestingly, the EM algorithm converges soon afterward 
as well, showing that for this problem EM spends little time 
in the region of parameter space in which a local analysis is 
valid. 



(b) 

Figure 2: Experimental results for the estimation of the 
parameters of a two-component Gaussian mixture (cf. 
Fig. 1). The separation of the Gaussians is half the 
separation in Fig. 1. 



(a) 



(b) 



Figure 3: The largest eigenvalues of the matrices 
H, E T HE, and E T PHE plotted as a function of the 
number of iterations. The plot in (a) is for the experi- 
ment in Fig. 1; (b) is for the experiment reported in Fig. 
2. 



the shape of / becomes more irregular. The condition 
number ^(0^)] increases to 352, n[E T H{Q^)E] in- 
creases to 216, and k[E t P (0^)11 (Q^E] increases to 
61. We see once again a significant improvement in the 
case of EM, by factors of 4 and 6. 

Fig. 3 shows that the matrix P has also reduced the 
largest eigenvalues of the Hessian from between 2000 to 
3000 to around 1. This demonstrates clearly the sta- 
ble convergence that is obtained via EM, without a line 
search or the need for external selection of a learning 
stepsize. 

In the remainder of the paper we provide some theo- 
retical analyses that attempt to shed some light on these 
empirical results. To illustrate the issues involved, con- 
sider a degenerate mixture problem in which the mixture 
has a single component. (In this case a.\ = 1.) Let us fur- 
thermore assume that the covariance matrix is fixed (i.e., 
only the mean vector m is to be estimated). The Hes- 
sian with respect to the mean m is H = —AS -1 and the 
EM projection matrix P is T,/N . For gradient ascent, we 
have k[E t HE] = k[S _1 ], which is larger than one when- 
ever £ ^ cl. EM, on the other hand, achieves a condi- 
tion number of one exactly (k[E t PHE] = k[PH] = 
k[I] = 1 and \ M [E T PHE] = 1). Thus, EM and New- 



ton's method are the same for this simple quadratic 
problem. For general non-quadratic optimization prob- 
lems, Newton retains the quadratic assumption, yield- 
ing fast convergence but possible divergence. EM is 
a more conservative algorithm that retains the conver- 
gence guarantee but also maintains quasi-Newton be- 
havior. We now analyze this behavior in more detail. 
We consider the special case of estimating the means in 
a Gaussian mixture when the Gaussians are well sepa- 
rated. 

Theorem 2 Consider the EM algorithm in Eq. (3), 
where the parameters oij and T,j are assumed to be 
known. Assume that the K Gaussian distributions are 
well separated, such that for sufficiently large k the pos- 
terior probabilities h- (t) are nearly zero or one. For 
such k, the condition number associated with EM is al- 
ways smaller than the condition number associated with 
gradient ascent. That is: 

K[E T p(e (k) )H(e (k) )E] < K [E T H(e (k) )E]. 

Furthermore, \m[E t P(Q^')H(Q^')E~\ approaches one 
as k goes to infinity. 



that the theorem can be extended to apply more widely, 
in particular to the case of the full EM update in which 
the mixing proportions and covariances are estimated, 
and also, within limits, to cases in which the means are 
not well separated. To obtain an initial indication as to 
possible conditions that can be usefully imposed on the 
separation of the mixture components, we have stud- 
ied the case in which the second term in Eq. (19) is 
neglected only for Ha and is retained for the Hij com- 
ponents, where j ^ i. Consider, for example, the case 
of a univariate mixture having two mixture components. 
For fixed mixing proportions and fixed covariances, the 
Hessian matrix (Eq. 17) becomes: 



H 



/in h 12 

h'21 h'22 



and the projection matrix (Eq. 19) becomes: 
P = 



-Ki 1 o 



where 



-h:, 



N 



JsyE^O. «■= 1.2 



2(k) 



Proof. The Hessian is 



H 



#n 

#21 



#12 
#22 



#2K 



Hri H_ 



K2 



H 



KK 



where 
#17 



d 2 i 



dm; dm,: 



(17) 



(18) 



N 



WT'EWW + Pn 



<*)\-i 



N 



Ti ( Uk),-i 



^(^(^-mj^-mif]^) 
t=i 

with jij(x^) = (S {j - hf\t))hf\t). The projection 
matrix P is 



^) = d 1 ag[^,.--^f ) 



KK)> 



where 



y(k) 
p (k) _ ±i_ h (k) () 

Z^t=i 

Given that h- (t)(l — h- it)) is negligible for suffi- 
ciently large k, the second term in Eq. (19) can be 
neglected, yielding Ha = -(SJ ) _1 J2t=i hj \t) and 
H = diag[#ii, • • • , Hkk]- This implies that PH = —I, 
and thus k[PH] = 1, whereas k[H] ^ 1. d 

This theorem, although restrictive in its assumptions, 
gives some indication as to why the projection matrix 
in the EM algorithm appears to condition the Hessian, 
yielding improved convergence. In fact, we conjecture 



1 N 

% = iP(i)E( 1 -' i f ) ( 1 ))''f ) w( i(i) -™)) T (^ i) -™.) i 



for i zfz j = 1, 2. If H is negative definite, (i.e., h\\\i22 — 
^12^21 < 0), then we can show that the conclusions of 
Theorem 2 remain true, even for Gaussians that are not 
necessarily well-separated. The proof is achieved via the 
following lemma: 

Lemma 1 Consider the positive definite matrix 



en 

<7"21 



Cl2 
C22 



For the diagonal matrix B = diag[<7 1:L , <t 22 ], we have 
k[5S] <«[£]. 



Proof. The eigenvalues of £ are the roots of (en 
A)(<7 2 2 - A) - (T 2 i(Ti2 = 0, which gives 



Am 



<7ll + (722 + 7 



(711 + (722 - 7 



7 = V( cr ll + cr 22) 2 - 4((7ii(7 2 2 - (721(7l 2 ) 



and 



«[£] 



<7ll + (722 + 7 
(711 + (722 - 7 



The condition number k[E] can be written as k[E] 
(1 + s)/(l — s) = /(s), where s is defined as follows: 



4((7ll(7 2 2 - (721(712) 



((711 + (722) 2 



Furthermore, the eigenvalues of BY, are the roots 

Of (1 - A)(l - A) - (<721<7i2)/(<Tii(T22) = 0, which 

gives A M = 1 + \/(o" 2iO"i2)/(o"iiO"22) and A m = 
1 - V( cr 2iO'i2)/(o'iiO'22)- Thus, defining r = 

\/(o-2io-i2)/(o-iia 2 2), we have k[5S] = (1 + r)/(l - r) = 
/(r). 

We now examine the quotient s/r: 



4(1 -r 2 ) 



((Tn +CT22) 2 /(o-llO-22) 

Given that (en + <722) 2 /(ciiC22) > 4, we have - > 

-\/l — (1 — r 2 ) = 1. That is, s > r. Since f(x) = 
(l + x)/(l — x) is a monotonically increasing function for 
x > 0, we have /(s) > /(r). Therefore, k[5S] < k[Y]. 
D 

We think that it should be possible to generalize 
this lemma beyond the univariate, two-component case, 
thereby weakening the conditions on separability in The- 
orem 2 in a more general setting. 

6 Conclusions 

In this paper we have provided a comparative analysis 
of algorithms for the learning of Gaussian mixtures. We 
have focused on the EM algorithm and have forged a link 
between EM and gradient methods via the projection 
matrix P. We have also analyzed the convergence of 
EM in terms of properties of the matrix P and the effect 
that P has on the likelihood surface. 

EM has a number of properties that make it a par- 
ticularly attractive algorithm for mixture models. It en- 
joys automatic satisfaction of probabilistic constraints, 
monotonic convergence without the need to set a learn- 
ing rate, and low computational overhead. Although EM 
has the reputation of being a slow algorithm, we feel 
that in the mixture setting the slowness of EM has been 
overstated. Although EM can indeed converge slowly 
for problems in which the mixture components are not 
well separated, the Hessian is poorly conditioned for 
such problems and thus other gradient-based algorithms 
(including Newton's method) are also likely to perform 
poorly. Moreover, if one's concern is convergence in like- 
lihood, then EM generally performs well even for these 
ill-conditioned problems. Indeed the algorithm provides 
a certain amount of safety in such cases, despite the poor 
conditioning. It is also important to emphasize that 
the case of poorly separated mixture components can 
be viewed as a problem in model selection (too many 
mixture components are being included in the model), 
and should be handled by regularization techniques. 

The fact that EM is a first order algorithm certainly 
implies that EM is no panacea, but does not imply that 
EM has no advantages over gradient ascent or superlin- 
ear methods. First, it is important to appreciate that 
convergence rate results are generally obtained for un- 
constrained optimization, and are not necessarily indica- 
tive of performance on constrained optimization prob- 
lems. Also, as we have demonstrated, there are condi- 
tions under which the condition number of the effective 



Hessian of the EM algorithm tends toward one, showing 
that EM can approximate a superlinear method. Finally, 
in cases of a poorly conditioned Hessian, superlinear con- 
vergence is not necessarily a virtue. In such cases many 
optimization schemes, including EM, essentially revert 
to gradient ascent. 

We feel that EM will continue to play an important 
role in the development of learning systems that empha- 
size the predictive aspect of data modeling. EM has in- 
deed played a critical role in the development of hidden 
Markov models (HMM's), an important example of pre- 
dictive data modeling. 4 EM generally converges rapidly 
in this setting. Similarly, in the case of hierarchical mix- 
tures of experts the empirical results on convergence in 
likelihood have been quite promising (Jordan & Jacobs, 
1994; Waterhouse k Robinson, 1994). Finally, EM can 
play an important conceptual role as an organizing prin- 
ciple in the design of learning algorithms. Its role in this 
case is to focus attention on the "missing variables" in 
the problem. This clarifies the structure of the algorithm 
and invites comparisons with statistical physics, where 
missing variables often provide a powerful analytic tool. 
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